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Trapped electron mode turbulence is studied by gyrokinetic simulations with the 
GYRO code and an analytical model including the effect of a poloidally varying elec- 
trostatic potential. Its impact on radial transport of high-Z trace impurities close to 
the core is thoroughly investigated and the dependence of the zero-flux impurity den- 
sity gradient (peaking factor) on local plasma parameters is presented. Parameters 
such as ion-to-electron temperature ratio, electron temperature gradient and main 
species density gradient mainly affect the impurity peaking through their impact on 
mode characteristics. The poloidal asymmetry, the safety factor and magnetic shear 
have the strongest effect on impurity peaking, and it is shown that under certain 
scenarios where trapped electron modes are dominant, core accumulation of high-Z 
impurities can be avoided. We demonstrate that accounting for the momentum con- 
servation property of the impurity-impurity collision operator can be important for 
an accurate evaluation of the impurity peaking factor. 



I. INTRODUCTION 

Turbulence driven by unstable drift waves is considered to be responsible for most of 
the observed cross-field particle and heat transport in the core of tokamaks. In particular, 
ion gyro-radius scale drift waves destabilized by the non-adiabatic response of trapped 
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electrons, the so-called trapped electron (TE) modes, can play an important role, specifi- 
cally in conditions where the electron heating power is large compared to the ion heating 
power and the electron temperature is larger than the ion temperature. They can also 
be important in transport barrier regions, where the density gradient is large [l|. 

Since its original discovery {2I, TE modes have been the topic of theoretical investi- 
gations. They are usually categorized into the dissipative and coUisionless classes 3| . The 
dissipative TE mode requires a strong temperature gradient and large collisionality, while 
the coUisionless TE mode - which is more likely to be destabilized in reactor relevant con- 
ditions - is driven by the electron curvature drift resonance and can be destabilized even 
in the absence of collisionality. The coUisionless TE mode can be driven purely by the 
main species density gradient, or by the electron temperature gradient. Consequently it 
is customary to further divide the coUisionless TE mode into density gradient driven and 
electron temperature gradient driven categories. The stability and the turbulent fluxes 



driven by TE modes have been analyzed in [4, l5|. It has been shown in that a 

quasilinear electrostatic approximation might retain much of the relevant physics of TE 
mode driven transport as it appears in nonlinear gyrokinetic simulations. The purpose of 
this paper is to study the impurity transport driven by TE modes. 

It is well known that accumulation of impurities - particularly those with high charge 
number - in the core of fusion plasmas has debilitating effect on fusion reactivity due to 
radiative losses and plasma dilution. Results of fluid and gyrokinetic simulations |7H23| 
indicated that the anomalous impurity transport driven by electrostatic microinstabilities 
in general, and TE modes in particular, is determined by the competition of three main 
mechanisms: curvature, thermo diffusion and parallel compressibility. The first of these 
contribute to an inward impurity transport (when the magnetic shear is positive), while 
thermodiffusion depends on the direction of the mode propagation, being inward for modes 
propagating in the electron diamagnetic direction such as the TE modes, however this 
contribution is negligible for high-Z impurities. The sign of the parallel compressibility 
contribution also depends on the direction of the mode propagation, but is instead outward 
for TE modes and has a charge to mass ratio dependence. 

In recent years attention has been directed towards the role of TE modes in impurity 
transport in plasmas with radio frequency (RF) heating. Various experiments reported 
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24|-|26| ■ In particular, it was shown 



reduced impurity accumulation in such circumstances 

in Ref. {2^] that impurity transport was more affected by the change in the plasmapa- 
rameters due to RF heating than by the generated sawtooth activity. In Ref. {3, Isj it 
was argued that in a TE-dominated ASDEX-U discharge with Electron Cyclotron Res- 
onance Heating (ECRH) the outward flows due to parallel compressibility explained the 
reduction in the impurity density peaking. In Ion Cyclotron Resonance Heating (ICRH) 
discharges on the JET tokamak the differences in the steady state impurity density profiles 
under minority heating (peaked impurity profiles) and mode conversion heating (hollow 
or flat impurity profiles) were partially explained by ITG and density gradient driven 
TE dominated transport, respectively [o], U8]. However, to be fully consistent with the 
observations an assumption of a sub-dominant electron temperature gradient driven TE 
mode was necessary. Further experimental studies in JET plasmas 26j showed a favor- 
able impact of ICRH in preventing the accumulation of metallic impurities accumulating 
in the core. However, in this case a theoretical explanation based on the presence of a 
TE mode driving outward impurity flux would be unsatisfactory since these plasmas were 
ITG dominated. 

Recently, a new possibility has emerged from the work reported in 29l-l3l| where the 
observed outward directed impurity flux is explained as an effect of poloidal asymmetries 
generated by the ICRH. The temperature anisotropy due to ICRH will trap the minority 
heated ions on the low field side, leading to the establishment of a poloidally varying 
"equilibrium" (i.e. non- fluctuating) electrostatic potential (such asymmetries have exper- 
imentally been demonstrated in |32|). The associated E x B drift acting as another degree 
of freedom for impurities to respond to electrostatic perturbations modifies the fluctuating 
impurity distribution. It has been shown that under experimentally relevant conditions 
the contribution of these E x B drifts to the impurity particle transport can be outward 
and might dominate the resulting steady state impurity gradient even in ITG dominant 
regimes. 

The density peaking of high-Z impurities in density or temperature gradient driven 
TE dominated plasmas under RF induced poloidal asymmetries is yet to be analyzed; this 
is the aim of the present paper. Apart from numerical simulations with gyro 33[ (mainly 
linear simulations, but a few nonlinear simulations are also performed for comparison) an 
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analytical model including the effect of poloidal asymmetries [3ll l3J] is utilized. The 
model is based on a solution to the linearized gyrokinetic equation, and it suggests that 
the impurity velocity pinch is governed by three separate contributions: one related to 
the magnetic drifts (combined effects of curvature and thermodiffusion pinch), another 
to the parallel impurity velocity pinch, and a third part arising due the E x B drift in 
a poloidally varying equilibrium electrostatic potential. Using this model, we present a 
systematic comparison of impurity transport driven by density and temperature driven 
trapped electron modes, highlighting the effect of the parallel impurity motion, collisions, 
magnetic geometry (shear and safety factor) and poloidal asymmetries. 

The remainder of the paper is organized as follows. In Sec. Ull we describe the 
baseline density and temperature gradient driven TE mode cases and the linear stability 
characteristics of them. In Sec. Illlt the density peaking of high-Z trace impurities is 
analyzed, and the dependence on relevant plasma parameters, such as electron density 
and temperature gradients, ion-to-electron temperature ratio, safety factor and magnetic 
shear, is presented. Impurity peaking factors are calculated in cases where the impurities 
are poloidally symmetrically distributed, but also in cases where a poloidally varying 
potential is present. The results are discussed and summarized in Sec. IIVI 



II. STABILITY 



The TE mode instability is driven by the electron logarithmic temperature gradient, 
a/LTei and/or the logarithmic density gradients, a/Ln-, whereas ITG modes are driven 
by the ion logarithmic temperature gradient a/LTi- Here Lna = — [dilnria) /dr]~'^ and 
L-Ta = — [c^(1iiTq,) /dr]~'^ represent the density and temperature scale lengths of particle 
species a, respectively, and a the outermost minor radius of the plasma. 

In this paper we will study two baseline coUisionless TE mode cases: one driven by 
the density gradients and one driven by the electron temperature gradient. For the second 
case the ion temperature gradient is set to zero, in order to obtain pure TE turbulence. 
This represents a situation with dominant central electron heating. Our baseline cases 
have the following local profile and magnetic geometry parameters: 

Case 1: Density gradient driven TE mode 
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Rq/q = 3, ro/a = 0.5, g = 2, s = 1, /3 = 0, a/L„ = 3, a/Lxe = o-lLxi = cl/Ltz = 1, 
Te = Ti = T^, z>ei = 0, p^o/a = 0.0035 

Case 2: Electron temperature gradient driven TE mode 

Ro/a = 3, ro/a = 0.375, q = 1.4, s = 0.8, (3 = 0, a/Ln = 1, a/Lre = 7/3, 
o/Ltz = 7/3, a/Ln = 0, Te = Tj = T,, z>ei = 0, pso/a = 0.0035 

Here the indices represent electrons (e), main ions (z) and impurities (z). The density 
gradient driven case (Case 1) is one of the GYRO standard cases in the GYRO non- 



linear gyrokinetic simulation database [35|. The electron temperature gradient driven 
case (Case 2) have been used in the fluid simulations presented in Ref. In both 

cases fully ionized nickel, Z = 28, is introduced in trace (i.e. Zuz/n^ <^ 1) quantities 
nz/ue = 2 X 10~^, however note that Z'^Uz/rie ~ C> {1) which is important for the ap- 
proximate model of the impurity peaking factor we will use. The use of nickel will ease 



comparison with previous work e.g. |9|, |3l|, but the main conclusions will be valid for any 
high-Z impurity. Ro is the major radius of the magnetic axis and tq the local reference 
minor radius, q is the safety factor and s = {r /q){dq/dr) the magnetic shear, while /3 
represents the ratio of plasma pressure to magnetic pressure. We note that electromag- 
netic fluctuations appearing for finite (3 have negligible effect on TE modes as trapped 
electrons cannot carry parallel current. In the main part of the paper, the plasma is 
assumed to be hot enough for collisions to be ignored and consequently the electron-ion 
collision frequency is Uei = 0, except when it is stated otherwise. 

This pap er considers turbulent fluxes. Neoclassical simulations of the baseline cases 
with NEO 36[ using u^i = 0.0058 c^/a (corresponding to Tg = 7 keV, = 3 x 10^^ m~'^. 
In A = 17 and a = 1 m) result in fluxes that are an order of magnitude smaller than the 
turbulent fluxes from nonlinear GYRO simulations. 

Nonlinear GYRO simulations of the baseline cases show that the largest fluxes occur 
in the vicinity of kgps = 0.15 for both of them; see Fig. [1] showing the poloidal wave number 
spectra of the gyro-Bohm normalized electron energy fluxes. Here kg is the poloidal wave- 
number and ps = Pso (1 + ecos6') the ion sound Larmor radius, where pso denotes ps at 
Ro, e = to/Rq is the inverse aspect ratio and 6 the extended poloidal angle. Consequently 
kgps = 0.15 is used in the quasilinear simulations for both cases. Note, however, that 
the maximum of the linear growth rates 7 are located at higher kgps as shown in Fig. |2l 
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As expected for a temperature gradient driven TE mode, the frequency increases, i.e. it 
propagates faster in the electron diamagnetic direction, with increasing kg in Case 2 {q]. 
Frequencies are given in Cg/a units, where = (Te/rrii)^^'^ is the ion sound speed and 
tq = a for the last closed flux surface. 

(TotoI.elec, energy) (TotoI.elec, energy) 
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FIG. 1: Normalized electron energy fluxes Qc/Qgb as functions of poloidal wave-number kgps 
from nonlinear GYRO simulations for CASE 1 (a) and CASE 2 (b). 
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FIG. 2: Linear growth rate 7 (circle markers, blue dashed lines) and real mode frequency oor 
(circle markers, red solid lines) as functions of poloidal wave-number kgps for CASE 1 (a) and 
Case 2 (b). Linear growth rate 7 (diamond markers, green dotted lines) and real mode frequency 
ujj. (diamond markers, orange dash-dotted lines) for the same cases but with parallel ion motion 
neglected in gyro. 



The perturbed electrostatic potential (p and eigenvalues u = + i'^ are obtained 
by linear electrostatic gyrokinetic initial- value calculations with gyro 33|. Linear initial- 
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value studies only consider the most unstable mode and any sub-dominant modes are 
neglected. In the simulations a model Grad-Shafranov magnetic equilibrium was used, 
where the O (e) corrections to the drift frequencies are retained. Flux-tube (periodic) 
boundary conditions were used, with a 128 point velocity space grid (8 energies, 8 pitch 
angles and two signs of velocity), the number of radial grid points is 6, and the number 
of poloidal grid points along particle orbits is 20 for trapped particles. The location of 
the highest energy grid point is at rriiv'^/ (2Tj) = 6. The ions were taken to be gyrokinetic 

1 /2 

and the electrons to be drift kinetic with the mass ratio (mj/me) = 60. 

The nonlinear electrostatic gyro simulations performed for the baseline cases also 
use gyrokinetic ions and drift kinetic electrons and the same velocity resolution as the 
linear simulations. At least 18 toroidal modes are used to model 1/4^^ of the torus, with 
the highest resolved poloidal wave number being kgps ~ 0.9. The number of radial grid 
points is 200. The simulations are run with the integration time step At = O.Ola/c^ for 
t > 200a/c,. 

Introducing a small collision rate is expected to have stabilizing effect on the colli- 
sionless TE mode, because the trapped electrons, driving the instability, can be detrapped. 
One of the most interesting distinctions between the two different branches of the TE 
modes we study concerns the dependence of the linear growth rate on coUisionality. If the 
TE mode is mainly driven by the electron temperature gradient, the mode is corapletely 
stabilized by collisions at a very low collision frequency, as was pointed out in Ref . j37| . As 
shown in Fig. [3l this was verified also in our simulations, where Case 2 was suppressed 
already for O^i > 0.015 c^/a while Case 1 persisted even for very high collisionalities, 
also consistent with earlier studies of density gradient driven TE modes Q (note 
that the z/gj-ranges plotted are different, and that Ur is positive for modes propagating to 
the electron diamagnetic direction according to GYRO conventions). For O^i > 1.0 Cg/a 
Case 1 even exhibits an increase in growth rate with increasing coUisionality. This could 
indicate that the mode is turning into a dissipative TE mode, but it can also suggest 
that TE modes driven by density gradients remain unstable even at large collisionalities. 
Further it can be noted that neglecting parallel ion motion in the simulations, leads to a 
small to moderate reduction of the real mode frequency while the growth rate is almost 
unaffected. 
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FIG. 3: Linear growth rate 7 (circle markers, blue dashed lines) and real mode frequency ujr 
(circle markers, red solid lines) as functions of electron-ion collision frequency Uei for Case 1 (a) 
(note the logarithmic i>ej-axis) and CASE 2 (b). Linear growth rate 7 (diamond markers, green 
dotted lines) and real mode frequency uJr (diamond markers, orange dash-dotted lines) for the 
same cases but with parallel ion motion neglected in gyro. 

Figure m shows how the perturbed potential varies with the extended poloidal angle. 
We see that both cases exhibit highly ballooned structures, concentrated to 6' G [— 7r,7r], 
and that there is no significant difference between including and not including parallel ion 
motion (note that the effect of parallel ion motion is expected to be stronger in cases of 
lower keps). 

III. IMPURITY DENSITY PEAKING 



for trace impurities 



34j . where the effect 



In this section, the zero flux density gradient (peaking factor _ 
is analyzed. We utilize a semi-analytical model introduced in Ref. 
of a poloidally varying equilibrium electrostatic potential is included. The focus is on 
the poloidally varying part of the electrostatic potential, and effects caused by a radial 
electric field, such as toroidal rotation, are neglected. However, we note, that the Coriolis 
drift or the poloidal redistribution of impurities due to centrifugal forces can have a 



non- negligible influence on impurity transport, as found in recent works |38|-|41|. The 
poloidally varying potential introduces an E x B drift frequency labeled w^; in the GK 
equation (Eq. ([3])), which disappears, ue = 0, in the symmetric case (note that gyro only 
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FIG. 4: Linear parallel mode structure of the perturbed potential (p {9) for CASE 1 (a) and 
Case 2 (b). Real part (red solid lines) and imaginary part (blue dashed lines) of (p. Real part 
(orange dash-dotted lines) and imaginary part (green dotted lines) of (/) for the same cases but 
with parallel ion motion neglected in gyro. Note that the actual resolution of the simulation 
covers O/tt = (—7,5), by gyro convention. 

considers the poloidally symmetric case). Poloidal variation can be caused by the presence 



of a species with strong temperature anisotropy [32| , which is the case in discharges with 
radio frequency (RF) heating of minority ions on the outboard side 42|, |43 |. 

The poloidally varying potential is assumed to be weak in the sense that eA^^j/T^ ^ 
1 (where is the temperature of species a). This implies that the effect of poloidal 
asymmetries on the main species can be neglected. This justifies the use of GYRO simula- 
tions neglecting poloidal asymmetries to obtain linear mode characteristics. By requiring 
Z ^ 1, we allow ZeA(j)E/Tz ~ O (1) and consequently the impurities can be poloidally 
asymmetrically distributed. Hence their E x B drift in the poloidally varying electrostatic 



potential 0s is not negligible. This model was presented in Ref. SjJ] and we refer to this 



work for more details. T^ 
given by Eq. (11) in Ref. 



le model we use for the equilibrium electrostatic potential is 



3l| 



Ze(j)E/Tz = —Kcos {9 — 6) , (1) 

where S represents the angular position where the impurity density has its maximum and 
K, sets the strength of the poloidal asymmetry. Thus the impurity density will be assumed 
to vary according to {6, r) = Uzo (?") A/'(^) with J^{0) = exp [k cos {6 — 6)] . In the model 
for ion cyclotron resonance heating (ICRH) driven asymmetries presented in Ref. j^sj 
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5 = TT is obtained, however, since impurity accumulation has also been observed at other 
poloidal locations we shall consider 5 = and 5 = 7r/2 cases as well. 



Note that, in contrast to Ref. 



3l|, the work presented here as well as in Ref. 34 1 



retains the effects of the parallel ion streaming in the GK equation (Eq. ([3])). In the 
case of the transport driven by TE mode turbulence, as we will show here, this term can 
significantly affect the impurity peaking. 



A. Zero flux impurity density gradient 



We consider particle transport driven by a single, representative, toroidal mode. The 
impurity peaking factor is calculated by requiring the linear impurity flux to vanish 



= (r.) 



B 



(2) 



where (■) denotes the flux surface average, [■] denotes imaginary part, is the per- 
turbed impurity density, gz the non-adiabatic part of the perturbed impurity distribution 
function, Jq is the Bessel function of the flrst kind, = kj_v±/uJcz, uJcz = ZeB/rriz is the 
cyclotron frequency and k± = {1 + s^O"^)^^"^ kg. Furthermore is the impurity mass, </>* 
is the complex conjugate of the perturbed electrostatic potential, and B is the strength 
of the equilibrium magnetic fleld. The subscripts y and ± denote the parallel and perpen- 
dicular directions with respect to the magnetic fleld. 

The non-adiabatic perturbed impurity distribution is obtained from the linearized 
GK equation 



^'11 dgz 



qR 89 



(^Dz - (^e)9z - C[gz 



£,li 



-i^^^{uj-ujJ^)(j)Joizz), (3) 



where a; = cu,. + ^7 is the mode frequency, /^q = n^ol^^z/SvrTj)^/^ exp(— f'/Tj.) is the 
equilibrium Maxwellian distribution function, E = rrizv'^ /2+Ze(j)E is the total unperturbed 
energy, fi = rrizV^/ {2B) is the magnetic moment, n^(r) = nzoexp[—Ze(f)E{r)/Tz] is the 
poloidally varying impurity density and Uzo is a flux function. The diamagnetic frequency 



is deflned as a;*2 = —keTz/ZeBLnz and = [1 



[X 



3/2) Lnz/Lrz], and x = v/vtz 



represents velocity normalized to the thermal speed vtz = {2Tz/mzY^'^ ■ The magnetic 
drift frequency is ujdz = —2k0Tz{x\/2+xVjV (9) / {nizUJczR), where V (9) = cos 9+s9 sin 9. 
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The E X B drift frequency of the particles in the equihbrium electrostatic field ue is 

kes6d(j)E 
B r o9 



and was derived in Appendix A of Ref. 311 (here, the d(j)E/dr part is dropped). C[-] is 
the collision operator. 

A solution to Eq. (|3]), and the subsequent expression for the peaking factor, is pre- 



sented in Ref. 3J]. It is a perturbative solution in the small parameter <^ 1, 

keeping terms up to 0{Z~^) in the expansion of g^. This is based on the fact that 
ujozjuj, uj^^juj, and Jo(^2) — 1 ~ are all ~ XjZ small, and that our ordering 

Ze(f)E/Tz ~ requires that ue/oj also is formally ~ 1/Z small. The solution as- 

sumes that impurity self-collisions dominate over collisions with unlike species, which 
follows from the ordering n^Z'^ jn^ ~ O (1), and the self-collisions are modeled by the full 
linearized impurity-impurity collision operator ci'i, maintaining the conservation prop- 
erties and self-adjointness. As earlier mentioned it is also assumed that </> and uj are 
known from the solution of the linear gyrokinetic-Maxwell system (obtained from GYRO) 
and that they are unaffected by the presence of trace impurities, and in particular their 
poloidal asymmetry. The impurity transit frequency v\\l (qR) is typically much smaller 
than the mode frequency u, and therefore magnetic (and electrostatic) trapping of the 
impurities can be neglected. 



The expression for the impurity peaking factor, a/Z/°^, is given in Eq. (8) of Ref. [34 1 
and can be modified into 



a CL a 2a^ Zrrij c, Ur 



2—{V)^ + -SK{esm{e-5)) 



de 



2 



/\<pn , (5) 



where (. . . )0 = (. . . A/'|0p)/(A/'|0p). To find this expression 0(e) corrections together 
with finite values of the mode eigenfunction outside the range [— tt, tt] of the extended 
poloidal angle have been neglected. As a consequence the expression is not valid in cases 
of highly elongated ballooning eigenfunctions, but as shown in Fig. H] the TE modes we 
study have a localized to this interval. It is interesting to note that up to the considered 
order, 0{Z~^), both finite Larmor radius (FLR) effects and the effects of collisions do 
not appear. Furthermore we see that cl/LP^z consists of three terms: the first term of 
Eq. ([5]) represents the contribution of the magnetic drift, the second term stems from the 
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E X B drift and is only non-zero when there is a poloidally varying potential, and the last 
term arises because of the irapurity parallel dynamics. The first two terms were present 
already in Eq. (14) of Ref. 3lj, but in that expression parallel ion/impurity dynamics 
was neglected. The term due to parallel dynamics contains only non-negative quantities, 
except ojr, and consequently impurity parallel dynamics acts to increase (decrease) the 
impurity peaking if u.^ is negative (positive). This leads to the conclusion that when 
the TE mode is the dominant instability {ur > 0) the parallel dynamics should act to 
decrease the impurity peaking, while the opposite is true for ITG modes. We note that 
in the poloidally symmetric (A/" = 1) case the parallel compressibility term in Eq. ([5]) 
is consistent with Eqs. (9) and (10) of Ref. j3] in the high Z limit when fcy is defined 
as (gi?o)~^(|f^e0P)/(|0P)- The effect of parallel dynamics was pointed out already in 
Ql, where it was observed that in conditions of strong electron temperature gradient the 
mechanism can be large enough to reverse the total pinch of trace impurities from inwards 
to outwards. 



B. Parametric dependences in poloidally symmetric cases 

In the present section we show how the impurity peaking factor and mode eigenval- 
ues depend on electron temperature gradient, ion-to-electron temperature ratio, electron 
density gradient and safety factor, when 0£; is poloidally symmetric. Results are pre- 
sented for both Case 1 and CASE 2, when parallel ion/impurity dynamics are included 
as well as when they are neglected. Results from Eq. (jS]) are compared to simulations by 
GYRO. Although the analytical model for the peaking factor given in Eq. ([5]) is based on 
a single linear mode, for completeness we will present a few cases where peaking factors 
have been determined from nonlinear simulations with gyro. The reason for this is to 
provide further insight into the validity of the analytical model. Note, that the analytical 
results from Eq. are not expected to agree exactly even with the linear results of 
GYRO since the approximation only retains effects up to order 1/Z as well as neglects 
0{e) corrections. Note that for cases with parallel ion/impurity dynamics neglected, the 
mode characteristics come from GYRO simulations where parallel compressibility effects 
have been turned off. Therefore they differ slightly from the magnetic drift contribution 
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in cases where parallel dynamics is included. This can be observed in e.g. Fig. \5^-h, 
comparing the blue dotted line with the orange dashed line. 



1. Temperature and temperature gradient dependences 

Figure |5^-b shows how the nickel peaking factor varies with electron temperature 
gradient and Fig. [St-d the corresponding eigenvalues. For both cases there is a slight 
increase in peaking for increasing a/ Lrpe-, but the dependence is generally very weak. This 
is what is expected from Eq. (ISl), where there is no explicit dependence on a/ Ltc- Instead 
the variations are caused by changes in mode frequency, shown in Fig. [5b-d, and perturbed 
potential through the parallel compressibility term, which is reflected in the fact that 
almost no variation at all is observed for the results with parallel compressibility neglected. 
Both the linear growth rate and real mode frequency increase with electron temperature 
gradient, where the increase in growth rate is expected since electron temperature gradient 
is one of the drives for TE modes. 

The peaking factor dependence on ion-to-electron temperature ratio (note that 
Tj/Te = Tz/Te) is illustrated in Fig. |6t-b. We see that it exhibits a similar dependence as 
to a/ Ltci which is weak and only enters through changes in uj^-, 7 and 0. How uo^ and 7 
are affected by changes in Ti/T^. is shown in Fig. [Gt-d. Case 1 shows a significant increase 
in growth rate with Ti/T^, and decrease in real mode frequency, while Case 2 is almost 
unaffected. This implies that, as Tj/Te increases. Case 1 is approaching an ITG mode. 
Similar trends were reported in terms of mode frequencies and growth rates in Ref. [iS^ 
where temperature and density gradient driven TE mode dominated plasmas were com- 
pared (having Ohmic and electron cyclotron heating, respectively); in the density gradient 
driven case Ur (7) was found to decrease (increase) with increasing Tj/Te, while weaker 
temperature ratio dependences of Ur and 7 were observed in the temperature gradient 
driven case. 

Earlier studies of TE modes have shown that the line ar g rowth rate is expected 



to decrease with increasing Tj/Tg, as seen in Fig. 5 of Ref. 4J], which is seemingly in 



disagreement with what is observed here. However, it has to be noted that the simulation 



in Ref. 



44| was performed keeping kgpi fixed, while here kgps is fixed. Since keps 
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FIG. 5: (a, b) Impurity peaking factor for trace nickel as function of electron temperature 
gradient a/LTe for CASE 1 (a) and CASE 2 (b). Red solid line is the peaking factor from 
Eq. dni), orange dashed line the magnetic drifts contribution, and green dash-dotted line the 
parallel compressibility contribution. Blue dotted line is the peaking factor from Eq. ([5]) without 
parallel compressibility effects. Red diamonds and blue dots correspond to gyro results with 
and without parallel compressibility effects respectively. 

(c, d) Linear growth rate 7 (circle markers, blue dashed lines) and real mode frequency ujr (circle 
markers, red solid lines) as functions of electron temperature gradient a/LTe for Case 1 (c) and 
Case 2 (d). Linear growth rate 7 (diamond markers, green dotted lines) and real mode frequency 
ujr (diamond markers, orange dash-dotted lines) for the same cases but with parallel ion motion 
neglected in gyro. 

(Te/Ti)^^"^ k^pi, a parametric scan over Ti/Tg keeping kgpi fixed, results in varying k^p^ 
accordingly and thus these scalings are not comparable. Furthermore in Ref. j4J] results 
are presented with 7 normalized to vxi/Ln oc while here we normalize to c^/a oc 
-y/T^. A test was performed, varying kgps to keep kgpi fixed, and a similar decrease with 
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FIG. 6: (a, b) Impurity peaking factor for trace nickel as function of ion-to-electron temperature 
ratio Ti/Te (note that Tj/Tg = T^/Tg) for Case 1 (a) and Case 2 (b). Red solid line is the 
peaking factor from Eq. ([5]), orange dashed line the magnetic drifts contribution, and green 
dash-dotted line the parallel compressibility contribution. Blue dotted line is the peaking factor 
from Eq. ([5]) without parallel compressibility effects. Red diamonds and blue dots correspond 
to GYRO results with and without parallel compressibility effects respectively, 
(c, d) Linear growth rate 7 (circle markers, blue dashed lines) and real mode frequency (circle 
markers, red solid lines) as functions of ion-to-electron temperature ratio T^/Tg for Case 1 (c) 
and Case 2 (d). Linear growth rate 7 (diamond markers, green dotted lines) and real mode 
frequency ujr (diamond markers, orange dash-dotted lines) for the same cases but with parallel 
ion motion neglected in gyro. 



increasing Tj/Tg was found. Also Ref. 45| finds an increase in linear growth rate with 
increasing Tj/Tg in agreement with our observations. 

A similar rather weak dependence of the impurity peaking on a/L^e and Tj/Tg was 
found in gyrokinetic simulations by j^, as long as the most unstable mode remained the 
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same. As observed in Case 1 here, for density gradient driven TE modes the peaking 
factor is typically positive. For temperature gradient driven TE modes however it can be 
negative which is also found in Case 2 where it is close to zero or even below, as seen in 
Figs. |5]d and [6)d. 

There is a small discrepancy between the values found by GYRO and the values 
calculated from Eq. For CASE 1 Eq. systematically overestimates the magnetic 
drift contribution, while for CASE 2 discrepancies arise mainly due to the parallel com- 
pressibility term. Still we see that the effect of the parallel impurity dynamics is rather 
well explained by the approximate solution. 

From Figs. [SHHt-d we see that neglecting the parallel ion motion in the GK equation 
(f||/ (qR) dgz/dO = in Eq. ([3])) only has a minor impact on the TE mode frequencies. 



2. Density gradient dependence 

The scaling of the impurity density peaking with the main species density peaking 
is shown in Fig. [7ti-b, and the corresponding eigenvalues in Fig. [7b-d. It is interesting to 
note how the impurity peaking shows a significant decrease in Case 1 for a/Lne = 1-0, 
compared to the other values, which is almost solely due to the change in the factor 
Ut/ [uI + 7^) in the parallel compressibility term of Eq. ([5]) (the variation in </> {9) with 
a/Lne is rather weak). The change in magnitude of this factor is clearly observed in 
Fig. [7b, where for a/Lne = 1-0, ooj. and 7 are comparable in size while for the other points 
7 is significantly larger. For density gradient driven TE modes a strong reduction in the 



peaking factor towards weaker density gradients has prevously been reported rn Ref, M, 

using both quasilinear and nonlinear gene |47| simulations. In Case 2 oor/ {oof. + 7^) 
experience a small reduction with increasing a/Lne, and this is reflected in the impurity 
peaking factor which is increased because of the smaller size of the parallel compressibility 
term. For a/Lne = 0.5 this case has a negative peaking factor. 

For a/Lne = 1.0 we could expect that Case 1 changes to become a temperature 
gradient driven TE mode, since then the density and temperature gradients are compara- 
ble {a/Lxe = 1-0 for Case 1). It seems that for temperature gradient driven TE modes 
the magnitude of the ratio of Uj. compared to 7 is usually larger than for density gradient 
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FIG. 7: (a, b) Impurity peaking factor for trace nickel as function of electron density gradient 
a/Lne for Case 1 (a) and Case 2 (b). Red solid line is the peaking factor from Eq. ([5]), 
orange dashed line the magnetic drifts contribution, and green dash-dotted line the parallel 
compressibility contribution. Blue dotted line is the peaking factor from Eq. ([5]) without parallel 
compressibility effects. Red diamonds and blue dots correspond to gyro results with and 
without parallel compressibility effects respectively. 

(c, d) Linear growth rate 7 (circle markers, blue dashed lines) and real mode frequency ujr 
(circle markers, red solid lines) as functions of electron density gradient a/Lne for CASE 1 (c) 
and Case 2 (d). Linear growth rate 7 (diamond markers, green dotted lines) and real mode 
frequency uJr (diamond markers, orange dash-dotted lines) for the same cases but with parallel 
ion motion neglected in gyro. 



driven TE modes. This has the consequence that the parallel dynamics becomes more 
important in reducing the impurity density peaking for temperature gradient driven TE 
modes. Further, we note that as could be expected from Eq. (|5]) the peaking factors are 
unaffected by changes in a/ Lm if parallel compressibility is neglected. For both cases the 
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linear growth rate and real mode frequency increase with increasing a/Lne- The eigen- 
values do not change much by neglecting parallel compressibility, except in Case 1 with 
a/Lne = 1.0 where an ITG mode is found instead of a TE mode. 



3. Safety factor dependence 

From the last term of Eq. (|5]), it is expected that the influence of impurity parallel 
dynamics on impurity peaking factor is strongly reduced with increasing safety factor. 
This is also what is observed in simulations, where the parallel compressibility contribution 
is significantly more negative with a small safety factor and as a consequence the peaking 
factor is more strongly reduced (see Fig. (H^-b). Since the safety factor typically is higher 
closer to the edge, it could be expected that in TE mode dominated plasmas the impurity 
parallel dynamics will reduce the peaking factor more and more, the closer to the core we 
look, while the opposite effect is expected in ITG dominated plasmas, as is confirmed by 



simulations presented in Ref. |2l|. However, as discussed in |5[, in a study of coUisionless 
TE modes, it should be noted that the ke leading to the largest fluxes is approximately 
inversely proportional to g, which is an effect we miss by keeping kg constant in our linear 
GYRO simulations. Since in the parallel compressibility contribution of Eq. (|5]) there is 
a factor 1/ (q'^ko), it is reduced to l/g if the mode leading to the largest fluxes should 
be considered. Furthermore it can also be noted that although the mode frequencies are 
relatively independent of q (as shown in Fig. Et-d), they are not independent of kgps (see 
Fig. [2]) but can on a very crude estimate be expected to vary linearly with kePs around the 
range of kgps we analyze. Because of the factor Ur/ (w^ + 7^) ^l^o found in the parallel 
compressibility contribution, this would imply that the dependence on q is completely 
canceled for the linear analysis of the mode leading to the largest fluxes. 

Figure IH^-b also includes impurity peaking factors determined from nonlinear 
GYRO simulations. These were calculated by linear interpolation of the impurity fluxes 
from two nonlinear GYRO simulations with different impurity density gradient. In these 
simulations, to keep an optimal /c^iPs-resolution around the peak part of the nonlinear en- 
ergy and particle flux spectra, the spacing between the simulated toroidal mode numbers 
are changed from case to case while the total number of toroidal modes are held fixed. In 
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FIG. 8: (a, b) Impurity peaking factor for trace nickel as function of safety factor q for CASE 1 (a) 
and Case 2 (b). Red solid line is the peaking factor from Eq. ([5]), orange dashed line the magnetic 
drifts contribution, and green dash-dotted line the parallel compressibility contribution. Blue 
dotted line is the peaking factor from Eq. ([5|) without parallel compressibility effects. Red 
diamonds and blue dots correspond to GYRO results with and without parallel compressibility 
effects respectively, while black hollow squares are results from nonlinear gyro runs, 
(c, d) Linear growth rate 7 (circle markers, blue dashed lines) and real mode frequency Ur (circle 
markers, red solid lines) as functions of safety factor q for Case 1 (c) and Case 2 (d). Linear 
growth rate 7 (diamond markers, green dotted lines) and real mode frequency ojr (diamond 
markers, orange dash-dotted lines) for the same cases but with parallel ion motion neglected in 

GYRO. 



the nonlinear simulation for the density gradient driven case we see a significantly weaker, 
but still existent, g-dependence than what the fixed- fcg linear simulation predicts. This 
may be understood from the above reasoning about the shift in the peak of the turbulent 
spectrum. On the other hand, in Case 2 the trend is found to be similar to the linear 
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predictions, and interestingly, the g-dependence is even stronger in the nonhnear case. 

The difference between including and not-including parallel ion dynamics for lower 
values of q is more pronounced in the electron temperature gradient driven CASE 2 than 
in the density gradient driven CASE 1 (a trend, consistent with the results of Ref. 13|). 
This is mostly due to that the factor Ur/ (w^ + 7^) Eq. ([5]) is larger for Case 2 than for 
Case 1 which can be seen from Fig. [Ht-d. Accordingly we could expect that for modes 
where this factor is relatively large, the effect of parallel impurity motion on the impurity 
peaking factor is strengthened. The differences in the dependence on q can only come 
through the g-dependence of u and cf) (6), since all the other parameters appearing in the 
third term of Eq. [5] are the same in the two cases for a fixed q. 



C. Poloidally asymmetric case 

In earlier studies with a poloidally varying potential present, magnetic shear has 
3een emphasized as one of the most important parameters affecting the impurity peaking 
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34| . This can be understood from its explicit appearance in the E x B drift term of 
Eq. (ISl). Consequently this section will focus on how the impurity peaking factor varies 
with magnetic shear, and results will be presented for the poloidally asymmetric cases with 
K, = 0.5 and 6 = (out-in asymmetry), 6 = 7i/2 (up-down asymmetry) and 6 = it (in-out 
asymmetry) in Eq. ([1]). We will omit the parametric dependence on other parameters in 
this section, because of the structure of the E x B drift term of Eq. (jS]) and refer to the 
results of Sec. IIIIBi The inclusion of the E x B drift term typically leads to the addition 
of a constant to the peaking factor, for other scalings than with s. This is because the 
only way this term can change in these scalings is through |</)p, and it is not varying by 
much. 

Figure |9^-b shows how the peaking factor depends on magnetic shear for Case 1 and 
Case 2 in both the symmetric and asymmetric cases. For Case 1 the symmetric peaking 
factor is mainly governed by the contribution from the magnetic drifts, which increases 
with s. The contribution from parallel compressibility is relatively small, and it is not 
affected much by a change in s. When introducing the asymmetry, the peaking factor 
changes significantly because of the E x B drift term which then is non-zero. In the in-out 
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FIG. 9: (a, b) Impurity peaking factor for trace nickel as function of magnetic shear s for 
Case 1 (a) and Case 2 (b). Red solid line is the peaking factor from Eq. ([5]) in the symmetric 
case, green dash-dotted line corresponds to out-in asymmetry, orange dashed line corresponds 
to up-down asymmetry, and black dotted line corresponds to in-out asymmetry. Red diamonds 
correspond to GYRO results. 

(c, d) Linear growth rate 7 (circle markers, blue dashed lines) and real mode frequency LOr (circle 
markers, red solid lines) as functions of magnetic shear s for Case 1 (c) and Case 2 (d). Linear 
growth rate 7 (diamond markers, green dotted lines) and real mode frequency Ur (diamond 
markers, orange dash-dotted lines) for the same cases but with parallel ion motion neglected in 

GYRO. 



asymmetric case there is a strong decrease of the impurity peaking, while on the contrary 
in the out-in asymmetric case there is a strong increase. The peaking factor remains 
almost unaffected by an up-down asymmetry. Since the E x B drift term becomes larger 
in magnitude with increasing s, the difference between the symmetric and asymmetric 
peaking factor is also increased with s. The reason why the E x B drift term leads to 
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a reduction (an increase) of the peaking factor for inboard (outboard) impurity accumu- 
lation is because that for 6 = n {6 = 0) the term (^sin(6' — 6))^ is negative (positive). 
Note however that if s < 0, this term changes sign and leads to an increase (decrease) 
for inboard (outboard) impurity accumulation, as shown in 3l|. The peaking factor in 
Case 2 shows a similar behavior to Case 1, with the main difference being that the 
contribution from the parallel compressibility term is significantly larger. Because of this 
term, the peaking factor can be negative even in the symmetric case. 

Furthermore in Fig. ^-d it can be noted that the linear eigenvalues are not strongly 
affected by changes in magnetic shear, although for s > 1.0 there is a small stabilizing 
effect with increasing shear. Ref. 4J] reported on a stabilizing effect which was weaker 
if the density gradients were large. This is consistent with what is observed here, since 
Case 1 has a larger density gradient than Case 2 and the stabilizing effect is weaker 
in Case 1. An increasing positive magnetic shear can stabilize a TE mode through 
FLR effects, but it can also drive the mode more unstable by increasing bad magnetic 
curvature. The dependence of linear growth rate on shear is consequently not trivial. The 
eigenvalues do not change much if the parallel ion dynamics is neglected. 



D. Collisions 

The model represented by Eq. ([5]) models impurity self-collisions by the full linearized 
impurity-impurity collision operator Cz} , and it is found that up to the considered order, 
0{Z~^), the effect of collisions does not appear explicitly. It thus predicts that the 
only way for collisions to affect the impurity peaking, is through their impact on the 
mode characteristics. This leads to a conclusion that in reality collisions should have a 
relatively weak influence on the impurity peaking factor. As many of the easily accessible 
gyrokinetic tools employ non-momentum-conserving model operators, it is interesting to 
see whether or not the form of the collision operator affects the above result. 

In Appendix A we present an alternative model which uses the Lorentz (or "pitch- 
angle scattering") operator to model impurity self-collisions, but is otherwise similar to 
the perturbative solution represented by Eq. ([5]). The most striking difference between the 
two models, is the appearance of the factor 1/ (1 + (x) /u) in the contribution related 
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to parallel dynamics, where u^, {x) [oc z/gj] is the deflection frequency. This implies that in 
the case of the Lorentz operator the effect of collisions appears explicitly in the expression 
for the impurity peaking factor, which is an artifact of C22[f||/2o] being different from 
zero. Figure [TO] shows a comparison between the two models for Case 1, but also for 
an ITG dominated case (earlier studied in 3J] also using fully ionized trace Nickel with 
local profile and geometry parameters: r/a = 0.3, Ro/a = 3, kgps = 0.3, q = 1.7, s = 1.5, 
a/Lne = 1-5, Tj/Tg = 0.85, a/Lxe = 2 and a/ L-n = 2.5). Case 2 is not considered since it 
is stabilized already at low coUisionality as shown in Fig. [31 The scalings illustrate that the 
impact of the Lorentz collisions starts to become important for z/gj > 0.1(~ |t<;|), and the 
two models start to diverge. The use of the Lorentz operator leads to an overestimation 
of the impurity peaking factor in the TE mode case, when the collision frequency is high, 
both in the poloidally symmetric case as well as in the asymmetric case. On the contrary, 
in the ITG mode case the use of the Lorentz operator leads to an underestimation of 
the peaking factor. This is expected because of the impact of 1/ (1 + iu^ {x) /u) in the 
parallel dynamics term, which decreased in magnitude with increasing u^, (x). Since 
parallel dynamics decreases (increases) the peaking factor for TE (ITG) modes we find 
an increase (a decrease) in the peaking factor with increasing ud {x). 



IV. CONCLUSIONS 



The paper presents a quasilinear study of two collisionless TE mode cases, driven by 
the density gradients, and the electron temperature gradient, respectively, including their 
mode characteristics and their effect on impurity transport. Mode characteristics has 
been obtained by linear gyrokinetic simulations using gyro. The poloidal wave number 
was chosen as kgps = 0.15 for both cases to represent the mode with the largest fluxes 
in nonlinear simulations. In agreement with previous studies, the electron temperature 
gradient driven mode is suppressed for small collisionalities, while the density gradient 
driven mode not only remains unstable even for very high collision frequencies, we also 
observe its transition to a dissipative TE mode. The dependence of growth rate and 
real frequency on safety factor and magnetic shear is non-monotonic and within small 
variations. 
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Uei for Case 1 (a) and an ITG dominated case (b) (note the logarithmic Pgi-axis). Red solid line is 
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from a model that utilize the Lorentz collision operator, and blue dashed line the corresponding 
in the in-out asymmetric case. 
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To investigate the peaking of high-Z trace impurities in tokamak plasmas we use an 
approximate gyrokinetic model and compare it to results obtained with gyro. It is ob- 
served that parameters such as Tj/Tg, a/Lxe and a/Lne mainly affect the peaking through 
their impact on mode characteristics, particularly the factor uor/ (w^ + 7^) is important to 
determine the effect of the impurity parallel dynamics. As noted before, in fluid modeling 
this factor enters directly into the parallel compressibility term of the approximate model, 
and is consequently responsible for determining the size of this contribution. Parameters 
describing the magnetic geometry, q and s, have a more significant influence on the peak- 
ing because of their explicit appearance in certain contributions. An increase in magnetic 
shear typically leads to an increase of the impurity peaking factor in the poloidally sym- 
metric case, because of the increase in the magnetic drift contribution. However in the 
poloidally asymmetric case since the term describing the E x B drift of impurities in the 
non-fluctuating electrostatic potential has an explicit linear shear dependence, an increase 
in shear can lead to a significant reduction or enhancement of the impurity peaking, de- 
pending on the location of the potential minimum. Increasing safety factor leads to a 
decrease of the relative significance of the impurity parallel dynamics contribution, but 
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the effect on the pealfing depends on the sign of and for TE modes, with Ur > 0, it 
results in an increase of the peaking factor. Nonhnear simulations in the density gradient 
driven TE case show only a very weak q scaling. This can be explained by a nonlinear 
shift in the poloidal wave number {kg ~ 1/q, as shown in [5]). However, in the tempera- 
ture gradient driven TE case, the q scaling is even stronger non-linearly than the fixed kg 
linear modeling predicts. 

The model using the conservation properties of the full linearized collision operator 
for impurity self collisions show that collisions can only indirectly affect impurity trans- 
port through changes in the mode characteristics. We show that when a non-momentum- 
conserving model operator is used, such as a Lorentz operator, the parallel compressibility 
contribution to the peaking factor is modified leading to errors in the coUisionality de- 
pendence. This effect becomes important when the impurity collision frequency becomes 
comparable to the mode frequency. 



Appendix A: Derivation peaking factor using a Lorentz collision operator 



In this section we derive a model for the impurity peaking factor similar to that 
represented by Eq. ([5]), but using the Lorentz collision operator instead of the linearized 
impurity-impurity collision operator. We neglect 0{e) corrections. 

The Lorentz collision operator for impurity self-collisions is given by 



I'D jx) d 

2 di 



(6) 



where is the deflection frequency for self-collisions ud (x) = z>22[Erf(x) — G{x)]/x^ 



v. 



\3/2 



and In A is the Coulomb logarithm. Erf (x) is the 
error-function and G (x) = [Erf (x) — xErf (x)] / (2a;^) the Chandrasekhar function. In 
the Lorentz operator ^ = x\\/x denotes the cosine of the pitch-angle. 

We assume the ordering udz/<^ ~ ^Izl^ ~ w^/w ~ Jq{zz) — 1 ~ and expand 
Qz in l/\fZ keeping terms up to 0{Z^^), i.e. gz ^ go + gi + g2- The O*'^ order solution of 
GK equation (Eq. [3]) is 

90 = — (7) 

z 

with C [(7o] = 0. This, added to the adiabatic response, —Ze(j)fzo/Tz, merely tells that 
the impurities are so heavy, and bound to the field lines through their high charge, that 
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they do not respond to electrostatic fluctuations to lowest order in \jZ. This justifles 
neglecting the effect of impurities on the mode characteristics, in spite of our assumption 
HzZ"^ /rif, ~ 1 (that is required to make self- collisions dominate). 
The 1*^* order GK equation reads 

n dgo 



iujgi - C[gi 



0, 



qR 89 

and is solved by assuming that the solution can be written in terms of Legendre poly- 
nomials P„ (0 as = g^x) Po (0 + 9l i^) Pi iO ■ Here Pq (0 = 1, (0 = ^ and 
we remind about the properties J^^d^PniO PmiO = 25m„/(2n + 1) and £[P„(^)] = 
—n {n + 1) Pn{C,) . The solution to Eq. ([8]) is found to be 

. v\i Zefzo d(j) 1 1 



qR Tz 86 u 1 + iuo/uj 
Note that for a momentum conserving collision operator C[gi oc i'n/^o] 
thus would not appear in Eq. (Q. 
The 2^^^^ order GK equation is 



(9) 



in Eq. 



iujg2 + C[g2] = iujozgo + 



n dgi Zef: 



zO 



I + 



(10) 



qR 89 T, 

The velocity anisotropies enter in udz and Zz-, which we can rewrite in terms of Legen- 
dre polynomials as udz = (1/3) [2Po (0 + P2 (0] ^d. and = (2/3) [Pq (0 - P2 (0] 4, 
where P2 (0 = (3^^ — 1) /2, and lodx = ^Dx{x), Zx = Zx{x) only depend on speed. 
Furthermore, by noting that 8£/86 = (where S = m^w^/2 + Ze(f)E) we can flnd the 
identity 



Ze8<i)i 



rriz 86 



■Po (0 . 



We now search for a solution to Eq. fIlOp of the form 

g2 = g'Ax)PoiO+9Ux)P2{0 



(12) 



where we realize that g"^ will not contribute to the particle flux since J^^d^P2 (0 — 0- 
By substituting Eq. (fT2l) into Eq. (fTOj) and collecting the parts proportional to Pq we flnd 
that 
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From g ^ Qo + gi + g2, only the g2 part of g2 contributes to the particle flux, thus the 
impurity peaking factor is found from solving 



= (r.) 



£vJo{z;,)g^(l)* 



3,, ^0 i* 



d'vg"2 



(14) 



where higher order than 1/Z corrections to the impurity flux are neglected. 
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